% 参数设置
t0 = 1;  % 初始时间
tf = 4;  % 终止时间
m = 18;  % 基函数的数量
N = 100; % 离散点的数量

% 时间变量变换
x0 = -1;
xf = 1;
c = (xf - x0) / (tf - t0);

% 生成离散点
x = linspace(x0, xf, N);
t = t0 + (tf - t0) * (x - x0) / (xf - x0);

% 初始化矩阵 A 和向量 b
A = zeros(N, m);
b = zeros(N, 1);

% 定义 Chebyshev 多项式基函数
for k = 1:m
    % 计算第 k 个 Chebyshev 多项式及其导数
    h = chebyshevT(k-1, x); % 第 k 个 Chebyshev 多项式
    h_prime = c * diff(chebyshevT(k-1, x)); % 一阶导数
    h_double_prime = c^2 * diff(diff(chebyshevT(k-1, x))); % 二阶导数
    
    % 确保导数数组大小与 t 一致
    h_prime = [h_prime, zeros(1,N - length(h_prime))]; % 补充零以匹配大小
    h_double_prime = [h_double_prime, zeros(1,N - length(h_double_prime))]; % 补充零以匹配大小
    
    % 填充矩阵 A
    A(:, k) = c^2 * t.^2 .* h_double_prime - c * t .* (t + 2) .* h_prime + (t + 2) .* h;
end

% 计算向量 b
f = zeros(N, 1); % 右侧函数 f(t) = 0
b = f;

% 应用初始条件
% y(1) = 1 和 y'(1) = 0
% 转换为 x = -1
g0 = 0; % g(x) 在 x = -1 处的值
g_prime_0 = 0; % g'(x) 在 x = -1 处的值

% 构建约束表达式
alpha = x .* (1 - x); % 选择 p(x) = 1 和 q(x) = x
beta = x.^2 / 2;

% 修改矩阵 A 和向量 b
A = A - alpha * g0 - beta * g_prime_0;
b = b - alpha * 1 - beta * 0; % y0 = 1, y'_0 = 0

% 求解最小二乘问题
coefficients = A \ b;

% 重建解
y = zeros(N, 1);
for k = 1:m
    y = y + coefficients(k) * chebyshevT(k-1, x);
end

% 添加约束项
y = y + alpha * 1 + beta * 0;

% 绘制结果
figure;
plot(t, y, 'LineWidth', 2);
xlabel('t');
ylabel('y(t)');
title('Least-Squares Solution of IVP');
legend('Least-Squares Solution');
grid on;